# PKO and foreign aid 

rm(list = ls())

library(foreign)
library(plm)
library(stargazer)
library(MatchIt)
library(cobalt)
library(car)

######################
# Part 1 analysis
######################

# import data
dt <- read.dta("pko-aid-main-data1.dta")

# Number of countries
length(unique(dt$gwno_loc))

# Descriptive statistics
summary(dt)
table(dt$pko_earlier)
sd(dt$ln_multi_aid)
sd(dt$pko_earlier)
sd(dt$ln_year_s_start)
sd(dt$intensity_level)
sd(dt$democracy_lag)
sd(dt$ln_cgdppc_lag)
sd(dt$ln_pop_lag)
sd(dt$ln_state_aid_lag)


# Table 1 model 1
t1mod1 <- plm(ln_multi_aid ~ pko_earlier + ln_multi_aid_lag, 
              data = dt, 
              index = c("gwno_loc","year"), model=c("within"))

stargazer(t1mod1, type = "text")

# Table 1 model 2
t1mod2 <- plm(ln_multi_aid ~ pko_earlier + ln_year_s_start +
                intensity_level + democracy_lag + 
                ln_cgdppc_lag + ln_pop_lag + ln_multi_aid_lag + 
                ln_state_aid_lag,
              data = dt, 
              index = c("gwno_loc","year"), model=c("within"))

stargazer(t1mod2, type = "text")

# VIF check
t1mod2.r <- lm(ln_multi_aid ~ pko_earlier + ln_multi_aid_lag + 
                 ln_state_aid_lag + ln_year_s_start + democracy_lag + 
                 ln_cgdppc_lag + ln_pop_lag + intensity_level,
               data = dt)
vif(t1mod2.r)
mean(vif(t1mod2.r))

# matching
colnames(dt)[colnames(dt)=="ln_year_s_start"] <- "Year since start(ln)"
colnames(dt)[colnames(dt)=="intensity_level"] <- "Intensity level"
colnames(dt)[colnames(dt)=="democracy_lag"] <- "Democracy"
colnames(dt)[colnames(dt)=="ln_cgdppc_lag"] <- "GDP pc(ln)"
colnames(dt)[colnames(dt)=="ln_pop_lag"] <- "Population(ln)"

m.out1 <- matchit(pko_earlier ~  
                    `Year since start(ln)` + `Intensity level` +
                    Democracy + 
                    `GDP pc(ln)` + `Population(ln)`, 
                  data = dt,
                  method = "cem")

bal.tab(m.out1,m.threshold=0.1)

m.data1 <- match.data(m.out1)
m.sum <- summary(m.out1)
m.sum


plot(m.sum, var.order = "unmatched", abs = FALSE, position = "topleft")


# Table 1 model 3
t1mod3 <- plm(ln_multi_aid ~ pko_earlier + ln_multi_aid_lag, 
              data = m.data1, 
              index = c("gwno_loc","year"), model=c("within"))

stargazer(t1mod3, type = "text")

#----------------
# Appendix Table C: Remove UN related multilateral aid
#----------------
rm(list = ls())
dt <- read.dta("pko-aid-main-data1.dta")
# Table C model 1
ACmod1 <- plm(ln_multi_aid_nonun ~ pko_earlier + ln_multi_aid_nonun_lag, 
              data = dt, 
              index = c("gwno_loc","year"), model=c("within"))

stargazer(ACmod1, type = "text")

# Table C model 2
ACmod2 <- plm(ln_multi_aid_nonun ~ pko_earlier + ln_year_s_start +
                intensity_level + democracy_lag + 
                ln_cgdppc_lag + ln_pop_lag + ln_multi_aid_nonun_lag + 
                ln_state_aid_lag,
              data = dt, 
              index = c("gwno_loc","year"), model=c("within"))

stargazer(ACmod2, type = "text")

#--------------------------------#
# Appendix D: Without lagged DV
#--------------------------------#

# Table D model 1
ADmod1 <- plm(ln_multi_aid ~ pko_earlier, 
              data = dt, 
              index = c("gwno_loc","year"), model=c("within"))

stargazer(ADmod1, type = "text")

# Table D model 2
ADmod2 <- plm(ln_multi_aid ~ pko_earlier + ln_year_s_start +
                intensity_level + democracy_lag + 
                ln_cgdppc_lag + ln_pop_lag +  
                ln_state_aid_lag,
              data = dt, 
              index = c("gwno_loc","year"), model=c("within"))

stargazer(ADmod2, type = "text")

#--------------------------------#
# Appendix E: Removing aid outflows that are controlled by the contributing donor governments
#--------------------------------#
rm(list = ls())
dt <- read.dta("pko-aid-main-data1.dta")

# Table E model 1
AEmod1 <- plm(ln_multi_aid2 ~ pko_earlier + ln_multi_aid2_lag, 
              data = dt, 
              index = c("gwno_loc","year"), model=c("within"))

stargazer(AEmod1, type = "text")

# Table E model 2
AEmod2 <- plm(ln_multi_aid2 ~ pko_earlier + ln_year_s_start +
                intensity_level + democracy_lag + 
                ln_cgdppc_lag + ln_pop_lag + ln_multi_aid2_lag + 
                ln_state_aid_lag,
              data = dt, 
              index = c("gwno_loc","year"), model=c("within"))

stargazer(AEmod2, type = "text")

#--------------------------------#
# Appendix F: Without fixed effects
#--------------------------------#

# Table F model 1
AFmod1 <- plm(ln_multi_aid ~ pko_earlier, 
              data = dt, 
              index = c("gwno_loc","year"), model=c("pooling"))

stargazer(AFmod1, type = "text")

# Table F model 2
AFmod2 <- plm(ln_multi_aid ~ pko_earlier + ln_year_s_start +
                intensity_level + democracy_lag + 
                ln_cgdppc_lag + ln_pop_lag +  
                ln_state_aid_lag,
              data = dt, 
              index = c("gwno_loc","year"), model=c("pooling"))

stargazer(AFmod2, type = "text")

#--------------------------------#
# Appendix G: Post conflict comparison
#--------------------------------#
rm(list = ls())
dt <- read.dta("pko-aid-main-data1.dta")
dt <- dt[dt$intensity_level==0, ]
# Table G model 1
AGmod1 <- plm(ln_multi_aid ~ pko_earlier + ln_multi_aid_lag, 
              data = dt, 
              index = c("gwno_loc","year"), model=c("within"))

stargazer(AGmod1, type = "text")

# Table G model 2
AGmod2 <- plm(ln_multi_aid ~ pko_earlier + ln_year_s_start +
                democracy_lag + 
                ln_cgdppc_lag + ln_pop_lag + ln_multi_aid_lag + 
                ln_state_aid_lag,
              data = dt, 
              index = c("gwno_loc","year"), model=c("within"))

stargazer(AGmod2, type = "text")

#--------------------------------#
# Appendix H: Post conflict comparison (snapshot ver)
#--------------------------------#
rm(list = ls())
dt <- read.dta("pko-aid-postconf.dta")
sub <- dt
sub <- unique(sub[sub$time_s_conf_term==1, c("gwno_loc","ln_multi_aid","conflictep_id")])
colnames(sub)[2] <- "ln_multi_aid1"
d5 <- dt[dt$time_s_conf_term==5, ]
newd <- merge(d5,sub,by = c("gwno_loc","conflictep_id"), all.x = T)
newd$ln_aid_diff <- newd$ln_multi_aid - newd$ln_multi_aid1
newd <- na.omit(newd)

# Table H model 1
AHmod1 <- lm(ln_aid_diff ~ pko_earlier, 
             data = newd)

stargazer(AHmod1, type = "text")

AHmod2 <- lm(ln_aid_diff ~ pko_earlier + 
                democracy_lag + 
                ln_cgdppc_lag + ln_pop_lag + ln_multi_aid_lag + 
                ln_state_aid_lag,
              data = newd)

stargazer(AHmod2, type = "text")


#######################
# Part 2 analysis 
#######################
rm(list = ls())
dt <- read.dta("pko-aid-main-data3.dta")

# Number of countries
length(unique(dt$gwno_loc))

# Descriptive statistics
summary(dt)
table(dt$govcap_elec_earlier)
sd(dt$v2clrspct)
sd(dt$govcap_elec_earlier)
sd(dt$ln_multi_aid_lag)
sd(dt$ln_year_s_start)
sd(dt$intensity_level)
sd(dt$democracy_lag)
sd(dt$ln_cgdppc_lag)
sd(dt$ln_pop_lag)
sd(dt$ln_state_aid_lag)


# Table 2 model 1
t2mod1 <- plm(v2clrspct ~ govcap_elec_earlier*ln_multi_aid_lag,
              data = dt, 
              index = c("gwno_loc","year"), model=c("within"))

stargazer(t2mod1, type = "text")


# Table 2 model 2
t2mod2 <- plm(v2clrspct ~ govcap_elec_earlier*ln_multi_aid_lag +
                ln_year_s_start + intensity_level + 
                democracy_lag + 
                ln_cgdppc_lag + ln_pop_lag + 
                ln_state_aid_lag,
              data = dt, 
              index = c("gwno_loc","year"), model=c("within"))

stargazer(t2mod2, type = "text")

# VIF check
t2mod2.r <- lm(v2clrspct ~ govcap_elec_earlier + ln_multi_aid_lag +
                 ln_state_aid_lag + 
                 ln_year_s_start + democracy_lag + 
                 ln_cgdppc_lag + ln_pop_lag + intensity_level,
               data = dt)
vif(t2mod2.r)
mean(vif(t2mod2.r))

# Appendix - include v2cltrnslw_lag
a2mod1 <- plm(v2clrspct ~ govcap_elec_earlier*ln_multi_aid_lag +
                ln_year_s_start + intensity_level + 
                democracy_lag + 
                ln_cgdppc_lag + ln_pop_lag + 
                ln_state_aid_lag + v2cltrnslw_lag,
              data = dt, 
              index = c("gwno_loc","year"), model=c("within"))

stargazer(a2mod1, type = "text")


# matching --------------
m.out2 <- matchit(govcap_elec_earlier ~  
                    ln_year_s_start + intensity_level + 
                    democracy_lag + 
                    ln_cgdppc_lag + ln_pop_lag + 
                    ln_state_aid_lag, 
                  data = dt,
                  method = "cem",
                  cutpoints = "q3")
m.data2 <- match.data(m.out2)
#write.dta(m.data2, "matched_part2.dta")

# Same thing with proper name
colnames(dt)[colnames(dt)=="ln_year_s_start"] <- "Year since start(ln)"
colnames(dt)[colnames(dt)=="intensity_level"] <- "Intensity level"
colnames(dt)[colnames(dt)=="democracy_lag"] <- "Democracy"
colnames(dt)[colnames(dt)=="ln_cgdppc_lag"] <- "GDP pc(ln)"
colnames(dt)[colnames(dt)=="ln_pop_lag"] <- "Population(ln)"
colnames(dt)[colnames(dt)=="ln_state_aid_lag"] <- "State aid(ln)"

m.out1 <- matchit(govcap_elec_earlier ~  
                    `Year since start(ln)` + `Intensity level` +
                    Democracy + `GDP pc(ln)` +
                    `Population(ln)` + `State aid(ln)`, 
                  data = dt,
                  method = "cem",
                  cutpoints = "q3")

bal.tab(m.out1,m.threshold=0.1)

m.data1 <- match.data(m.out1)
m.sum <- summary(m.out1)
# Figure 1
plot(m.sum, var.order = "unmatched", abs = FALSE, position = "topleft")

t2mod3 <- lm(v2clrspct ~ govcap_elec_earlier*ln_multi_aid_lag, 
             data = m.data1, weights = weights)

library(lmtest)
library(sandwich)

# Table 2 model 3
coeftest(t2mod3, vcov. = vcovHC)
stargazer(t2mod3, type = "text")

# Marginal effect (Results were obtained from Stata)
# Figure 2
rm(list = ls())

d <- read.table(text = '
 Capa/Elec tasks (0),   .0300953  
 Capa/Elec tasks (1),   .0503351', sep = ',')
d

par(mar=c(3,7,6,6), font=2, font.lab=2, cex=1)
plot(d$V2, xaxt = 'n', pch = 19, ylim = c(0,0.1),
     xlab = "", ylab = "The difference in the 
     predictive margins",
     main = "Marginal effects of multilateral aid (ln) 
     with 95% confidence intervals")
axis(1, at = seq(nrow(d)), labels = d$V1)
abline(h=c(0), col=c("black"), lty=c(2))
lines(c(1,1),c( .0025993  ,  .0575912))
lines(c(2,2),c( .0303231  , .0703472))


# Marginal effect (Results were obtained from Stata) from IV regression
# Figure 3
rm(list = ls())

d <- read.table(text = '
 Capa/Elec tasks (0),    -.0047999  
 Capa/Elec tasks (1),   .0633153', sep = ',')
d

par(mar=c(3,7,6,6), font=2, font.lab=2, cex=1)
plot(d$V2, xaxt = 'n', pch = 19, ylim = c(-0.02,0.1),
     xlab = "", ylab = "The difference in the 
     predictive margins",
     main = "Marginal effects of multilateral aid (ln) 
     with 95% confidence intervals")
axis(1, at = seq(nrow(d)), labels = d$V1)
abline(h=c(0), col=c("black"), lty=c(2))
lines(c(1,1),c(  -.015759  , .0061592))
lines(c(2,2),c( .0289063  ,  .0977244))


# Marginal effect (Results were obtained from Stata) from Sample selection
# Figure 4
rm(list = ls())

d <- read.table(text = '
 Capa/Elec tasks (0),    -.0242179  
 Capa/Elec tasks (1),   .0670663', sep = ',')
d

par(mar=c(3,7,6,6), font=2, font.lab=2, cex=1)
plot(d$V2, xaxt = 'n', pch = 19, ylim = c(-0.08,0.1),
     xlab = "", ylab = "The difference in the 
     predictive margins",
     main = "Marginal effects of multilateral aid (ln) 
     with 95% confidence intervals")
axis(1, at = seq(nrow(d)), labels = d$V1)
abline(h=c(0), col=c("black"), lty=c(2))
lines(c(1,1),c(   -.077883  ,  .0294472))
lines(c(2,2),c( .0377279  ,  .0964048))

#-------------------------#
# Appendix Table I (With additional control variable)
#-------------------------#
rm(list = ls())
dt <- read.dta("pko-aid-main-data3.dta")

# Table I model 1
AImod1 <- plm(v2clrspct ~ govcap_elec_earlier*ln_multi_aid_lag +
                v2cltrnslw_lag,
              data = dt, 
              index = c("gwno_loc","year"), model=c("within"))

stargazer(AImod1, type = "text")

# Table I model 2
AImod2 <- plm(v2clrspct ~ govcap_elec_earlier*ln_multi_aid_lag +
                ln_year_s_start + intensity_level + 
                democracy_lag + 
                ln_cgdppc_lag + ln_pop_lag + 
                ln_state_aid_lag + v2cltrnslw_lag,
              data = dt, 
              index = c("gwno_loc","year"), model=c("within"))

stargazer(AImod2, type = "text")

#-------------------------#
# Appendix Figure B (With additional control variable, IV regression)
#----------------------
rm(list = ls())
# Results were obtained from Stata
d <- read.table(text = '
 Capa/Elec tasks (0),    -.0033006   
 Capa/Elec tasks (1),    .0325131 ', sep = ',')
d

par(mar=c(3,7,6,6), font=2, font.lab=2, cex=1)
plot(d$V2, xaxt = 'n', pch = 19, ylim = c(-0.02,0.08),
     xlab = "", ylab = "The difference in the 
     predictive margins",
     main = "Marginal effects of multilateral aid (ln) 
     with 95% confidence intervals")
axis(1, at = seq(nrow(d)), labels = d$V1)
abline(h=c(0), col=c("black"), lty=c(2))
lines(c(1,1),c(   -.0128819  ,  .0062806))
lines(c(2,2),c(  .0014855  ,  .0635406))

#-------------------------#
# Appendix Figure C (With additional control variable, sample selection model)
#----------------------
rm(list = ls())
# Results were obtained from Stata
d <- read.table(text = '
 Capa/Elec tasks (0),     -.0216761    
 Capa/Elec tasks (1),    .018777  ', sep = ',')
d

par(mar=c(3,7,6,6), font=2, font.lab=2, cex=1)
plot(d$V2, xaxt = 'n', pch = 19, ylim = c(-0.04,0.04),
     xlab = "", ylab = "The difference in the 
     predictive margins",
     main = "Marginal effects of multilateral aid (ln) 
     with 95% confidence intervals")
axis(1, at = seq(nrow(d)), labels = d$V1)
abline(h=c(0), col=c("black"), lty=c(2))
lines(c(1,1),c(   -.0376037  , -.0057485))
lines(c(2,2),c(  .0012133  ,  .0363407))

#-------------------------#
# Appendix Table L (Without UN aid)
#-------------------------#
rm(list = ls())
dt <- read.dta("pko-aid-main-data3.dta")

# Table L model 1
ALmod1 <- plm(v2clrspct ~ govcap_elec_earlier*ln_multi_aid_nonun_lag,
              data = dt, 
              index = c("gwno_loc","year"), model=c("within"))

stargazer(ALmod1, type = "text")

# Table I model 2
ALmod2 <- plm(v2clrspct ~ govcap_elec_earlier*ln_multi_aid_nonun_lag +
                ln_year_s_start + intensity_level + 
                democracy_lag + 
                ln_cgdppc_lag + ln_pop_lag + 
                ln_state_aid_lag,
              data = dt, 
              index = c("gwno_loc","year"), model=c("within"))

stargazer(ALmod2, type = "text")

#-------------------------#
# Appendix Figure D (Without UN aid, IV regression)
#----------------------
rm(list = ls())
# Results were obtained from Stata
d <- read.table(text = '
 Capa/Elec tasks (0),     -.0127126      
 Capa/Elec tasks (1),    .0663458  ', sep = ',')
d

par(mar=c(3,7,6,6), font=2, font.lab=2, cex=1)
plot(d$V2, xaxt = 'n', pch = 19, ylim = c(-0.03, 0.11),
     xlab = "", ylab = "The difference in the 
     predictive margins",
     main = "Marginal effects of multilateral aid non-UN (ln) 
     with 95% confidence intervals")
axis(1, at = seq(nrow(d)), labels = d$V1)
abline(h=c(0), col=c("black"), lty=c(2))
lines(c(1,1),c(    -.0249183 ,  -.0005069))
lines(c(2,2),c(  .030644  ,  .1020475))

#-------------------------#
# Appendix Figure E (Without UN aid, Sample selection model)
#----------------------
rm(list = ls())
# Results were obtained from Stata
d <- read.table(text = '
 Capa/Elec tasks (0),     -.0001714        
 Capa/Elec tasks (1),   .0636352    ', sep = ',')
d

par(mar=c(3,7,6,6), font=2, font.lab=2, cex=1)
plot(d$V2, xaxt = 'n', pch = 19, ylim = c(-0.05, 0.1),
     xlab = "", ylab = "The difference in the 
     predictive margins",
     main = "Marginal effects of multilateral aid non-UN (ln) 
     with 95% confidence intervals")
axis(1, at = seq(nrow(d)), labels = d$V1)
abline(h=c(0), col=c("black"), lty=c(2))
lines(c(1,1),c(   -.0442457  ,   .043903))
lines(c(2,2),c( .0456609  ,  .0816096))

#-------------------------#
# Appendix Table O (Using multilateral aid 2)
#-------------------------#
rm(list = ls())
dt <- read.dta("pko-aid-main-data3.dta")

# Table O model 1
AOmod1 <- plm(v2clrspct ~ govcap_elec_earlier*ln_multi_aid2_lag,
              data = dt, 
              index = c("gwno_loc","year"), model=c("within"))

stargazer(AOmod1, type = "text")

# Table I model 2
AOmod2 <- plm(v2clrspct ~ govcap_elec_earlier*ln_multi_aid2_lag +
                ln_year_s_start + intensity_level + 
                democracy_lag + 
                ln_cgdppc_lag + ln_pop_lag + 
                ln_state_aid_lag,
              data = dt, 
              index = c("gwno_loc","year"), model=c("within"))

stargazer(AOmod2, type = "text")

#-------------------------#
# Appendix Figure F (Using multilateral aid2, IV regression)
#----------------------
rm(list = ls())
# Results were obtained from Stata
d <- read.table(text = '
 Capa/Elec tasks (0),     -.0068989        
 Capa/Elec tasks (1),   .0659664    ', sep = ',')
d

par(mar=c(3,7,6,6), font=2, font.lab=2, cex=1)
plot(d$V2, xaxt = 'n', pch = 19, ylim = c(-0.02, 0.11),
     xlab = "", ylab = "The difference in the 
     predictive margins",
     main = "Marginal effects of multilateral aid 2 (ln) 
     with 95% confidence intervals")
axis(1, at = seq(nrow(d)), labels = d$V1)
abline(h=c(0), col=c("black"), lty=c(2))
lines(c(1,1),c(   -.0171284  ,  .0033307))
lines(c(2,2),c(  .0312899 ,   .1006429))

#-------------------------#
# Appendix Figure G (Using multilateral aid2, Sample selection model)
#----------------------
rm(list = ls())
# Results were obtained from Stata
d <- read.table(text = '
 Capa/Elec tasks (0),    -.0226067         
 Capa/Elec tasks (1),   .0672047   ', sep = ',')
d

par(mar=c(3,7,6,6), font=2, font.lab=2, cex=1)
plot(d$V2, xaxt = 'n', pch = 19, ylim = c(-0.08, 0.1),
     xlab = "", ylab = "The difference in the 
     predictive margins",
     main = "Marginal effects of multilateral aid 2 (ln) 
     with 95% confidence intervals")
axis(1, at = seq(nrow(d)), labels = d$V1)
abline(h=c(0), col=c("black"), lty=c(2))
lines(c(1,1),c(   -.0777687  ,  .0325553))
lines(c(2,2),c( .0383081  ,  .0961013))

